2A: 1st Order Spatial Point Patterns Analysis

In this exercise, we will learn to analyze spatial point patterns in R, including importing geospatial data, performing kernel density estimation and nearest neighbor analysis, and visualizing results using spatstat, sf, and tmap packages.

Author

Teng Kok Wai (Walter)

Published

August 28, 2024

Modified

September 1, 2024

1 Exercise 2A Reference

R for Geospatial Data Science and Analytics - 4  1st Order Spatial Point Patterns Analysis Methods

2 Overview

Spatial Point Pattern Analysis is the evaluation of the pattern or distribution, of a set of points on a surface. The point can be location of:

  • events such as crime, traffic accident and disease onset, or
  • business services (coffee and fastfood outlets) or facilities such as childcare and eldercare.

Using appropriate functions of spatstat, this hands-on exercise aims to discover the spatial point processes of childecare centres in Singapore.

The specific questions we would like to answer are as follows:

  • are the childcare centres in Singapore randomly distributed throughout the country?
  • if the answer is not, then the next logical question is where are the locations with higher concentration of childcare centres?

3 Learning Outcome

  • Importing and managing geospatial data using sf and spatstat packages.
  • Converting simple feature data frames to various spatial data formats.
  • Performing kernel density estimation (KDE) using different methods and bandwidths.
  • Conducting nearest neighbor analysis to test spatial point distribution.
  • Visualizing spatial patterns and KDE results using tmap and spatstat.
  • Handling duplicate spatial points and creating point patterns for analysis.
  • Applying appropriate statistical methods for confirmatory spatial point pattern analysis.

Here’s the improved section with relevant information presented in a table format:

4 The Data

Dataset Description Source Format
CHILDCARE Point data containing location and attributes of childcare centers. Data.gov.sg GeoJSON
MP14_SUBZONE_WEB_PL Polygon data with URA 2014 Master Plan Planning Subzone boundaries. Data.gov.sg ESRI Shapefile
CoastalOutline Polygon data representing Singapore’s national boundary. Singapore Land Authority (SLA) ESRI Shapefile

This table summarizes the datasets used for the analysis, including their descriptions, sources, and formats.

Here’s the improved section with relevant information presented in a table format:

5 Installing and Loading the R Packages

The following R packages will be used in this exercise:

Package Purpose Use Case in Exercise
sf Imports, manages, and processes vector-based geospatial data. Handling vector geospatial data in R.
spatstat Provides tools for point pattern analysis. Performing 1st- and 2nd-order spatial point pattern analysis and deriving kernel density estimation (KDE).
raster Reads, writes, and manipulates gridded spatial data (raster). Converting image outputs from spatstat into raster format.
maptools Offers tools for manipulating geographic data. Converting spatial objects into ppp format for use with spatstat.
tmap Creates static and interactive thematic maps using cartographic quality elements. Plotting static and interactive point pattern maps.

To install and load these packages in R, use the following code:

pacman::p_load(sf, raster, spatstat, tmap, tidyverse)

6 Spatial Data Wrangling

6.1 Importing the Spatial Data

To import the three geographical datasets, we will use st_read() from sf.

childcare_sf <- st_read("data/child-care-services-geojson.geojson")
Reading layer `child-care-services-geojson' from data source 
  `/Users/walter/code/isss626/isss626-gaa/Hands-on_Ex/Hands-on_Ex02/data/child-care-services-geojson.geojson' 
  using driver `GeoJSON'
Simple feature collection with 1545 features and 2 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 103.6824 ymin: 1.248403 xmax: 103.9897 ymax: 1.462134
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
sg_sf <- st_read(dsn = "data", layer="CostalOutline")
Reading layer `CostalOutline' from data source 
  `/Users/walter/code/isss626/isss626-gaa/Hands-on_Ex/Hands-on_Ex02/data' 
  using driver `ESRI Shapefile'
Simple feature collection with 60 features and 4 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 2663.926 ymin: 16357.98 xmax: 56047.79 ymax: 50244.03
Projected CRS: SVY21
mpsz_sf <- st_read(dsn = "data",
                layer = "MP14_SUBZONE_WEB_PL")
Reading layer `MP14_SUBZONE_WEB_PL' from data source 
  `/Users/walter/code/isss626/isss626-gaa/Hands-on_Ex/Hands-on_Ex02/data' 
  using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21

6.2 Inspect and Reproject to Same Projection System

Before we can use these data for analysis, it is important for us to ensure that they are projected in same projection system.

First, we check the childcare dataset.

st_crs(childcare_sf)
Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

This dataset is using the WGS84 crs. We will reproject all the dataset to SVY21 crs for standardization and analysis.

childcare_sf <- st_transform(childcare_sf , crs = 3414)
st_crs(childcare_sf)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]

The childcare dataset has been reprojected to SVY21 successfully.

Next, we inspect the Coastal Outline dataset.

st_crs(sg_sf)
Coordinate Reference System:
  User input: SVY21 
  wkt:
PROJCRS["SVY21",
    BASEGEOGCRS["SVY21[WGS84]",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

Notice that this dataset is using SVY21 crs, however the ID provided is EPSG:9001 does not match the intended ID, EPSG:3414 of SVY21. In this case, we will set the crs to the correct ID using the code block below.

sg_sf <- st_set_crs(sg_sf, 3414)
st_crs(sg_sf)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]

Similarly, we will inspect the Master Plan Subzone Dataset.

st_crs(mpsz_sf)
Coordinate Reference System:
  User input: SVY21 
  wkt:
PROJCRS["SVY21",
    BASEGEOGCRS["SVY21[WGS84]",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

Since the ID is also EPSG:9001, we will set the crs to EPSG:3414 too.

mpsz_sf <- st_set_crs(mpsz_sf, 3414)
st_crs(mpsz_sf)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]

7 Mapping the Geospatial Datasets

After checking the referencing system of each geospatial data data frame, it is also useful for us to plot a map to show their spatial patterns.

7.1 Static Map

# add polygon layer of the coastal outline of sg island
tm_shape(sg_sf)+ tm_polygons() +
# add polygon layer of the subzone based on sg masterplan
tm_shape(mpsz_sf) + tm_polygons() +
# add dot layer to show the locations of childcare centres
tm_shape(childcare_sf) + tm_dots() +
tm_layout()

When all the 3 datasets are overlayed together, it shows the locations of childcare centres on the Singapore island. Since all the geospatial layers are within the same map context, it means their referencing system and coordinate values are referred to similar spatial context. This consistency is crucial for accurate geospatial analysis.

7.2 Interactive Map

Alternatively, we can also prepare a pin map by using the code block below.

tmap_mode('view')

# tm_basemap("Esri.WorldGrayCanvas") +
# tm_basemap("OpenStreetMap") +
tm_basemap("Esri.WorldTopoMap") +
tm_shape(childcare_sf) +
  tm_dots(alpha = 0.5)
tmap_mode('plot')

In interactive mode, tmap uses the Leaflet for R API, allowing you to freely navigate, zoom, and click on features for detailed information. You can also change the map’s background using layers like ESRI.WorldGrayCanvas, OpenStreetMap, and ESRI.WorldTopoMap, with ESRI.WorldGrayCanvas as the default.

Tip

Remember to switch back to plot mode after interacting to avoid connection issues and limit interactive maps to fewer than 10 in RMarkdown documents for Netlify publishing.

8 Geospatial Data Wrangling

While simple feature data frames are becoming more popular compared to sp’s Spatial* classes, many geospatial analysis packages still require data in the Spatial* format. This section will show you how to convert a simple feature data frame to sp’s Spatial* class.

8.1 Converting sf data frames to sp’s Spatial* class

The code block below uses as_Spatial() of sf package to convert the three geospatial data from simple feature data frame to sp’s Spatial* class.

childcare <- as_Spatial(childcare_sf)
mpsz <- as_Spatial(mpsz_sf)
sg <- as_Spatial(sg_sf)

After the sf dataframe to sp Spatial* conversion, let’s inspect the Spatial* classes.

childcare
class       : SpatialPointsDataFrame 
features    : 1545 
extent      : 11203.01, 45404.24, 25667.6, 49300.88  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 2
names       :    Name,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           Description 
min values  :   kml_1, <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>ADDRESSBLOCKHOUSENUMBER</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSBUILDINGNAME</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSPOSTALCODE</th> <td>018989</td> </tr><tr bgcolor=""> <th>ADDRESSSTREETNAME</th> <td>1, MARINA BOULEVARD, #B1 - 01, ONE MARINA BOULEVARD, SINGAPORE 018989</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSTYPE</th> <td></td> </tr><tr bgcolor=""> <th>DESCRIPTION</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>HYPERLINK</th> <td></td> </tr><tr bgcolor=""> <th>LANDXADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor="#E3E3F3"> <th>LANDYADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor=""> <th>NAME</th> <td>THE LITTLE SKOOL-HOUSE INTERNATIONAL PTE. LTD.</td> </tr><tr bgcolor="#E3E3F3"> <th>PHOTOURL</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSFLOORNUMBER</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>08F73931F4A691F4</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20200826094036</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSUNITNUMBER</th> <td></td> </tr></table></center> 
max values  : kml_999,                  <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>ADDRESSBLOCKHOUSENUMBER</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSBUILDINGNAME</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSPOSTALCODE</th> <td>829646</td> </tr><tr bgcolor=""> <th>ADDRESSSTREETNAME</th> <td>200, PONGGOL SEVENTEENTH AVENUE, SINGAPORE 829646</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSTYPE</th> <td></td> </tr><tr bgcolor=""> <th>DESCRIPTION</th> <td>Child Care Services</td> </tr><tr bgcolor="#E3E3F3"> <th>HYPERLINK</th> <td></td> </tr><tr bgcolor=""> <th>LANDXADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor="#E3E3F3"> <th>LANDYADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor=""> <th>NAME</th> <td>RAFFLES KIDZ @ PUNGGOL PTE LTD</td> </tr><tr bgcolor="#E3E3F3"> <th>PHOTOURL</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSFLOORNUMBER</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>379D017BF244B0FA</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20200826094036</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSUNITNUMBER</th> <td></td> </tr></table></center> 
mpsz
class       : SpatialPolygonsDataFrame 
features    : 323 
extent      : 2667.538, 56396.44, 15748.72, 50256.33  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 15
names       : OBJECTID, SUBZONE_NO, SUBZONE_N, SUBZONE_C, CA_IND, PLN_AREA_N, PLN_AREA_C,       REGION_N, REGION_C,          INC_CRC, FMEL_UPD_D,     X_ADDR,     Y_ADDR,    SHAPE_Leng,    SHAPE_Area 
min values  :        1,          1, ADMIRALTY,    AMSZ01,      N, ANG MO KIO,         AM, CENTRAL REGION,       CR, 00F5E30B5C9B7AD8,      16409,  5092.8949,  19579.069, 871.554887798, 39437.9352703 
max values  :      323,         17,    YUNNAN,    YSSZ09,      Y,     YISHUN,         YS,    WEST REGION,       WR, FFCCF172717C2EAF,      16409, 50424.7923, 49552.7904, 68083.9364708,  69748298.792 
sg
class       : SpatialPolygonsDataFrame 
features    : 60 
extent      : 2663.926, 56047.79, 16357.98, 50244.03  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 4
names       : GDO_GID, MSLINK, MAPID,              COSTAL_NAM 
min values  :       1,      1,     0,             ISLAND LINK 
max values  :      60,     67,     0, SINGAPORE - MAIN ISLAND 

The geospatial data have been converted into their respective sp’s Spatial* classes.

8.2 Converting the Spatial* Class into Generic sp Format

spatstat requires the analytical data in ppp object form. There is no straightforward way to convert a Spatial* classes into ppp object. We need to convert the Spatial classes* into Spatial object first.

Tip

ppp refers to planar point pattern. It is used to represent spatial point patterns in spatstat; it contains the event locations with possibly associated marks, and the observation window where the events occur.

see Chapter 18 The spatstat package | Spatial Statistics for Data Science: Theory and Practice with R

The codes block below converts the Spatial* classes into generic sp objects.

childcare_sp <- as(childcare, "SpatialPoints")
sg_sp <- as(sg, "SpatialPolygons")

Next, we can display the sp objects properties.

childcare_sp
class       : SpatialPoints 
features    : 1545 
extent      : 11203.01, 45404.24, 25667.6, 49300.88  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
sg_sp
class       : SpatialPolygons 
features    : 60 
extent      : 2663.926, 56047.79, 16357.98, 50244.03  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 

Let’s further inspect the differences between Spatial* classes and generic sp object with the example of childcare and childcare_sp object.

head(childcare)
   Name
1 kml_1
2 kml_2
3 kml_3
4 kml_4
5 kml_5
6 kml_6
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     Description
1 <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>ADDRESSBLOCKHOUSENUMBER</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSBUILDINGNAME</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSPOSTALCODE</th> <td>760742</td> </tr><tr bgcolor=""> <th>ADDRESSSTREETNAME</th> <td>742, YISHUN AVENUE 5, #01 - 470, SINGAPORE 760742</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSTYPE</th> <td></td> </tr><tr bgcolor=""> <th>DESCRIPTION</th> <td>Child Care Services</td> </tr><tr bgcolor="#E3E3F3"> <th>HYPERLINK</th> <td></td> </tr><tr bgcolor=""> <th>LANDXADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor="#E3E3F3"> <th>LANDYADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor=""> <th>NAME</th> <td>AVERBEL CHILD DEVELOPMENT CENTRE PTE LTD</td> </tr><tr bgcolor="#E3E3F3"> <th>PHOTOURL</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSFLOORNUMBER</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>AEA27114446235CE</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20200826094036</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSUNITNUMBER</th> <td></td> </tr></table></center>
2                                    <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>ADDRESSBLOCKHOUSENUMBER</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSBUILDINGNAME</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSPOSTALCODE</th> <td>159053</td> </tr><tr bgcolor=""> <th>ADDRESSSTREETNAME</th> <td>20, LENGKOK BAHRU, #02 - 05, SINGAPORE 159053</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSTYPE</th> <td></td> </tr><tr bgcolor=""> <th>DESCRIPTION</th> <td>Child Care Services</td> </tr><tr bgcolor="#E3E3F3"> <th>HYPERLINK</th> <td></td> </tr><tr bgcolor=""> <th>LANDXADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor="#E3E3F3"> <th>LANDYADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor=""> <th>NAME</th> <td>AWWA LTD.</td> </tr><tr bgcolor="#E3E3F3"> <th>PHOTOURL</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSFLOORNUMBER</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>86B24416FB1663C6</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20200826094036</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSUNITNUMBER</th> <td></td> </tr></table></center>
3        <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>ADDRESSBLOCKHOUSENUMBER</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSBUILDINGNAME</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSPOSTALCODE</th> <td>556912</td> </tr><tr bgcolor=""> <th>ADDRESSSTREETNAME</th> <td>22, LI HWAN VIEW, GOLDEN HILL ESTATE, SINGAPORE 556912</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSTYPE</th> <td></td> </tr><tr bgcolor=""> <th>DESCRIPTION</th> <td>Child Care Services</td> </tr><tr bgcolor="#E3E3F3"> <th>HYPERLINK</th> <td></td> </tr><tr bgcolor=""> <th>LANDXADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor="#E3E3F3"> <th>LANDYADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor=""> <th>NAME</th> <td>BABIES BY-THE-PARK PTE. LTD.</td> </tr><tr bgcolor="#E3E3F3"> <th>PHOTOURL</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSFLOORNUMBER</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>F971CBBA973E1AE5</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20200826094036</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSUNITNUMBER</th> <td></td> </tr></table></center>
4 <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>ADDRESSBLOCKHOUSENUMBER</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSBUILDINGNAME</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSPOSTALCODE</th> <td>569139</td> </tr><tr bgcolor=""> <th>ADDRESSSTREETNAME</th> <td>3, ANG MO KIO STREET 62, #01 - 36, LINK@AMK, SINGAPORE 569139</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSTYPE</th> <td></td> </tr><tr bgcolor=""> <th>DESCRIPTION</th> <td>Child Care Services</td> </tr><tr bgcolor="#E3E3F3"> <th>HYPERLINK</th> <td></td> </tr><tr bgcolor=""> <th>LANDXADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor="#E3E3F3"> <th>LANDYADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor=""> <th>NAME</th> <td>Baby Elk Infant Care Pte Ltd</td> </tr><tr bgcolor="#E3E3F3"> <th>PHOTOURL</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSFLOORNUMBER</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>86A4F25D1C7C9D85</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20200826094036</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSUNITNUMBER</th> <td></td> </tr></table></center>
5                           <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>ADDRESSBLOCKHOUSENUMBER</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSBUILDINGNAME</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSPOSTALCODE</th> <td>467961</td> </tr><tr bgcolor=""> <th>ADDRESSSTREETNAME</th> <td>22A, KEW DRIVE, SINGAPORE 467961</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSTYPE</th> <td></td> </tr><tr bgcolor=""> <th>DESCRIPTION</th> <td>Child Care Services</td> </tr><tr bgcolor="#E3E3F3"> <th>HYPERLINK</th> <td></td> </tr><tr bgcolor=""> <th>LANDXADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor="#E3E3F3"> <th>LANDYADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor=""> <th>NAME</th> <td>BABYPLANET MONTESSORI PTE. LTD.</td> </tr><tr bgcolor="#E3E3F3"> <th>PHOTOURL</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSFLOORNUMBER</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>CFE3F056F8171C7B</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20200826094036</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSUNITNUMBER</th> <td></td> </tr></table></center>
6                       <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>ADDRESSBLOCKHOUSENUMBER</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSBUILDINGNAME</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSPOSTALCODE</th> <td>598523</td> </tr><tr bgcolor=""> <th>ADDRESSSTREETNAME</th> <td>3 Jalan Kakatua, JURONG PARK, SINGAPORE 598523</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSTYPE</th> <td></td> </tr><tr bgcolor=""> <th>DESCRIPTION</th> <td>Child Care Services</td> </tr><tr bgcolor="#E3E3F3"> <th>HYPERLINK</th> <td></td> </tr><tr bgcolor=""> <th>LANDXADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor="#E3E3F3"> <th>LANDYADDRESSPOINT</th> <td>0</td> </tr><tr bgcolor=""> <th>NAME</th> <td>BAMBINI CHILDCARE LLP</td> </tr><tr bgcolor="#E3E3F3"> <th>PHOTOURL</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSFLOORNUMBER</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>2B4F0B285ED28C4A</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20200826094036</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSUNITNUMBER</th> <td></td> </tr></table></center>
head(childcare_sp)
class       : SpatialPoints 
features    : 1 
extent      : 27976.73, 27976.73, 45716.7, 45716.7  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 

Note that the Spatial* classes contain more attribute data as compared its generic sp object counterpart.

Tip

Differences between Spatial* classes and generic sp object

  • Data Storage: SpatialPoints stores only the coordinates, while SpatialPointsDataFrame stores both coordinates and additional attribute data

  • Functionality: SpatialPointsDataFrame allows for more complex operations and analyses due to the additional data it holds

see Introduction to spatial points in R - Michael T. Hallworth, Ph.D.

8.3 Converting the Generic sp Format into spatstat’s ppp Format

Now, we will use as.ppp() function of spatstat to convert the spatial data into spatstat’s ppp object format.

childcare_ppp <- as.ppp(childcare_sf)
childcare_ppp
Marked planar point pattern: 1545 points
marks are of storage type  'character'
window: rectangle = [11203.01, 45404.24] x [25667.6, 49300.88] units

Let’s examine the difference by plotting childcare_ppp:

plot(childcare_ppp)

We can also view the summary statistics of the newly created ppp object by using the code block below.

summary(childcare_ppp)
Marked planar point pattern:  1545 points
Average intensity 1.91145e-06 points per square unit

Coordinates are given to 11 decimal places

marks are of type 'character'
Summary:
   Length     Class      Mode 
     1545 character character 

Window: rectangle = [11203.01, 45404.24] x [25667.6, 49300.88] units
                    (34200 x 23630 units)
Window area = 808287000 square units
Tip

Be aware of the warning message regarding duplicates. In spatial point pattern analysis, duplicates can be a significant issue. The statistical methods used for analyzing spatial point patterns often assume that the points are distinct and non-coincident.

8.4 Handling Duplicated Points

We can check the duplication in a ppp object by using the duplicated function with different configurations.

Tip

Tips on using duplicated

If rule=“spatstat” (the default), two points are deemed identical if their coordinates are equal according to ==, and their marks are equal according to ==. This is the most stringent possible test. If rule=“unmark”, duplicated points are determined by testing equality of their coordinates only, using ==. If rule=“deldir”, duplicated points are determined by testing equality of their coordinates only, using the function duplicatedxy in the package deldir, which currently uses duplicated.data.frame. Setting rule=“deldir” will ensure consistency with functions in the deldir package.

see R: Determine Duplicated Points in a Spatial Point Pattern

# duplicated(childcare_ppp)
# any(duplicated(childcare_ppp))
rules <- c("spatstat", "deldir", "unmark")

duplicate_counts <- list()
for (rule in rules) {
  duplicates <- duplicated(childcare_ppp, rule = rule)
  num_duplicates <- sum(duplicates)
  duplicate_counts[[rule]] <- num_duplicates
}

print(duplicate_counts)
$spatstat
[1] 0

$deldir
[1] 0

$unmark
[1] 74

Using unmark, we can find 74 duplicates.

We can also use anyDuplicated(x) as it is a faster version of any(duplicated(x)).

sum(anyDuplicated(childcare_ppp, rule="unmark"))
[1] 0

To count the number of coincident points, we will use the multiplicity() function as shown in the code block below. see R: Multiplicity for more info.

multiplicity(childcare_ppp)

If we want to know how many locations have more than one point event:

sum(multiplicity(childcare_ppp) > 1)
[1] 0
# double check
coincident_points <- duplicated(childcare_ppp,  rule="deldir")
coincident_coordinates <- childcare_ppp[coincident_points]
print(coincident_coordinates)
Marked planar point pattern: 0 points
marks are of storage type  'character'
window: rectangle = [11203.01, 45404.24] x [25667.6, 49300.88] units

The output shows that there are 128(???) 0 duplicated point events.

tmap_mode('view')
tm_basemap("Esri.WorldTopoMap") +
tm_shape(childcare) +
  tm_dots(alpha=0.4,
          size=0.05)
tmap_mode('plot')

8.5 How to Spot Duplicate Points on the Map

There are three ways to overcome this problem. The easiest way is to delete the duplicates. But, that will also mean that some useful point events will be lost.

The second solution is use jittering, which will add a small perturbation to the duplicate points so that they do not occupy the exact same space.

The third solution is to make each point “unique” and then attach the duplicates of the points to the patterns as marks, as attributes of the points. Then you would need analytical techniques that take into account these marks.

8.5.1 Jittering

The code block below implements the jittering approach.

childcare_ppp_jit <- rjitter(childcare_ppp,
                             retry=TRUE,
                             nsim=1,
                             drop=TRUE)
any(duplicated(childcare_ppp_jit))
[1] FALSE

8.6 Creating owin Object

When analysing spatial point patterns, it is a good practice to confine the analysis with a geographical area like Singapore boundary. In spatstat, an object called owin is specially designed to represent this polygonal region.

The code block below is used to covert sg SpatialPolygon object into owin object of spatstat.

sg_owin <- as.owin(sg_sf)

The output object can be displayed by using plot() function:

plot(sg_owin)

And using summary() function of Base R:

summary(sg_owin)
Window: polygonal boundary
50 separate polygons (1 hole)
                 vertices         area relative.area
polygon 1 (hole)       30     -7081.18     -9.76e-06
polygon 2              55     82537.90      1.14e-04
polygon 3              90    415092.00      5.72e-04
polygon 4              49     16698.60      2.30e-05
polygon 5              38     24249.20      3.34e-05
polygon 6             976  23344700.00      3.22e-02
polygon 7             721   1927950.00      2.66e-03
polygon 8            1992   9992170.00      1.38e-02
polygon 9             330   1118960.00      1.54e-03
polygon 10            175    925904.00      1.28e-03
polygon 11            115    928394.00      1.28e-03
polygon 12             24      6352.39      8.76e-06
polygon 13            190    202489.00      2.79e-04
polygon 14             37     10170.50      1.40e-05
polygon 15             25     16622.70      2.29e-05
polygon 16             10      2145.07      2.96e-06
polygon 17             66     16184.10      2.23e-05
polygon 18           5195 636837000.00      8.78e-01
polygon 19             76    312332.00      4.31e-04
polygon 20            627  31891300.00      4.40e-02
polygon 21             20     32842.00      4.53e-05
polygon 22             42     55831.70      7.70e-05
polygon 23             67   1313540.00      1.81e-03
polygon 24            734   4690930.00      6.47e-03
polygon 25             16      3194.60      4.40e-06
polygon 26             15      4872.96      6.72e-06
polygon 27             15      4464.20      6.15e-06
polygon 28             14      5466.74      7.54e-06
polygon 29             37      5261.94      7.25e-06
polygon 30            111    662927.00      9.14e-04
polygon 31             69     56313.40      7.76e-05
polygon 32            143    145139.00      2.00e-04
polygon 33            397   2488210.00      3.43e-03
polygon 34             90    115991.00      1.60e-04
polygon 35             98     62682.90      8.64e-05
polygon 36            165    338736.00      4.67e-04
polygon 37            130     94046.50      1.30e-04
polygon 38             93    430642.00      5.94e-04
polygon 39             16      2010.46      2.77e-06
polygon 40            415   3253840.00      4.49e-03
polygon 41             30     10838.20      1.49e-05
polygon 42             53     34400.30      4.74e-05
polygon 43             26      8347.58      1.15e-05
polygon 44             74     58223.40      8.03e-05
polygon 45            327   2169210.00      2.99e-03
polygon 46            177    467446.00      6.44e-04
polygon 47             46    699702.00      9.65e-04
polygon 48              6     16841.00      2.32e-05
polygon 49             13     70087.30      9.66e-05
polygon 50              4      9459.63      1.30e-05
enclosing rectangle: [2663.93, 56047.79] x [16357.98, 50244.03] units
                     (53380 x 33890 units)
Window area = 725376000 square units
Fraction of frame area: 0.401

9 Combining Point Events Object and Owin Object

For the last step of geospatial data wrangling, we will extract childcare events that are located within Singapore by using the code block below.

childcareSG_ppp = childcare_ppp[sg_owin]

The output object combined both the point and polygon feature in one ppp object class as shown below.

summary(childcareSG_ppp)
Marked planar point pattern:  1545 points
Average intensity 2.129929e-06 points per square unit

Coordinates are given to 11 decimal places

marks are of type 'character'
Summary:
   Length     Class      Mode 
     1545 character character 

Window: polygonal boundary
50 separate polygons (1 hole)
                 vertices         area relative.area
polygon 1 (hole)       30     -7081.18     -9.76e-06
polygon 2              55     82537.90      1.14e-04
polygon 3              90    415092.00      5.72e-04
polygon 4              49     16698.60      2.30e-05
polygon 5              38     24249.20      3.34e-05
polygon 6             976  23344700.00      3.22e-02
polygon 7             721   1927950.00      2.66e-03
polygon 8            1992   9992170.00      1.38e-02
polygon 9             330   1118960.00      1.54e-03
polygon 10            175    925904.00      1.28e-03
polygon 11            115    928394.00      1.28e-03
polygon 12             24      6352.39      8.76e-06
polygon 13            190    202489.00      2.79e-04
polygon 14             37     10170.50      1.40e-05
polygon 15             25     16622.70      2.29e-05
polygon 16             10      2145.07      2.96e-06
polygon 17             66     16184.10      2.23e-05
polygon 18           5195 636837000.00      8.78e-01
polygon 19             76    312332.00      4.31e-04
polygon 20            627  31891300.00      4.40e-02
polygon 21             20     32842.00      4.53e-05
polygon 22             42     55831.70      7.70e-05
polygon 23             67   1313540.00      1.81e-03
polygon 24            734   4690930.00      6.47e-03
polygon 25             16      3194.60      4.40e-06
polygon 26             15      4872.96      6.72e-06
polygon 27             15      4464.20      6.15e-06
polygon 28             14      5466.74      7.54e-06
polygon 29             37      5261.94      7.25e-06
polygon 30            111    662927.00      9.14e-04
polygon 31             69     56313.40      7.76e-05
polygon 32            143    145139.00      2.00e-04
polygon 33            397   2488210.00      3.43e-03
polygon 34             90    115991.00      1.60e-04
polygon 35             98     62682.90      8.64e-05
polygon 36            165    338736.00      4.67e-04
polygon 37            130     94046.50      1.30e-04
polygon 38             93    430642.00      5.94e-04
polygon 39             16      2010.46      2.77e-06
polygon 40            415   3253840.00      4.49e-03
polygon 41             30     10838.20      1.49e-05
polygon 42             53     34400.30      4.74e-05
polygon 43             26      8347.58      1.15e-05
polygon 44             74     58223.40      8.03e-05
polygon 45            327   2169210.00      2.99e-03
polygon 46            177    467446.00      6.44e-04
polygon 47             46    699702.00      9.65e-04
polygon 48              6     16841.00      2.32e-05
polygon 49             13     70087.30      9.66e-05
polygon 50              4      9459.63      1.30e-05
enclosing rectangle: [2663.93, 56047.79] x [16357.98, 50244.03] units
                     (53380 x 33890 units)
Window area = 725376000 square units
Fraction of frame area: 0.401
plot(childcareSG_ppp)

10 First-order Spatial Point Patterns Analysis

In this section, you will learn how to perform first-order SPPA by using spatstat package. The hands-on exercise will focus on:

  • deriving kernel density estimation (KDE) layer for visualising and exploring the intensity of point processes,
  • performing Confirmatory Spatial Point Patterns Analysis by using Nearest Neighbour statistics.

10.1 Kernel Density Estimation

In this section, you will learn how to compute the kernel density estimation (KDE) of childcare services in Singapore.

10.1.1 Computing Kernel Density Estimation Using Automatic Bandwidth Selection Method

  • bw.diggle() automatic bandwidth selection method. Other recommended methods are bw.CvL(), bw.scott() or bw.ppl().
  • The smoothing kernel used is gaussian, which is the default. Other smoothing methods are: “epanechnikov”, “quartic” or “disc”.
  • The intensity estimate is corrected for edge effect bias by using method described by Jones (1993) and Diggle (2010, equation 18.9). The default is FALSE.
kde_childcareSG_bw <- density(childcareSG_ppp,
                              sigma=bw.diggle,
                              edge=TRUE,
                            kernel="gaussian")

The plot() function of Base R is then used to display the kernel density derived.

plot(kde_childcareSG_bw)

The density values of the output range from 0 to 0.000035 which is way too small to comprehend. This is because the default unit of measurement of SVY21 is in meter. As a result, the density values computed is in “number of points per square meter”.

Before we move on to next section, it is good to know that you can retrieve the bandwidth used to compute the kde layer by using the code block below.

bw <- bw.diggle(childcareSG_ppp)
bw
   sigma 
298.4095 

10.1.2 Rescaling KDE values

rescale.ppp() is used below to convert the unit of measurement from meter to kilometer:

childcareSG_ppp.km <- rescale.ppp(childcareSG_ppp, 1000, "km")

Now, we can re-run density() using the resale data set and plot the output kde map.

kde_childcareSG.bw <- density(childcareSG_ppp.km, sigma=bw.diggle, edge=TRUE, kernel="gaussian")
plot(kde_childcareSG.bw)

Since we just did a rescaling operation, the output image looks identical to the earlier version with the only changes in terms of data values.

10.2 Working with Different Automatic Bandwidth Methods

Beside bw.diggle(), there are 3 other spatstat functions can be used to determine the bandwidth, they are: bw.CvL(), bw.scott(), and bw.ppl().

Let us take a look at the bandwidth return by these automatic bandwidth calculation methods by using:

 bw.CvL(childcareSG_ppp.km)
   sigma 
4.543278 
bw.scott(childcareSG_ppp.km)
 sigma.x  sigma.y 
2.224898 1.450966 
bw.ppl(childcareSG_ppp.km)
    sigma 
0.3897114 
bw.diggle(childcareSG_ppp.km)
    sigma 
0.2984095 

Baddeley et al. (2016) suggest using the bw.ppl() algorithm, as it tends to produce more appropriate values when the pattern consists predominantly of tight clusters. However, they also note that if the aim of a study is to detect a single tight cluster amidst random noise, the bw.diggle() method is likely to be more effective.

To compare the output of using bw.diggle and bw.ppl methods:

kde_childcareSG.ppl <- density(childcareSG_ppp.km,
                               sigma=bw.ppl,
                               edge=TRUE,
                               kernel="gaussian")
par(mfrow=c(1,2))

plot(kde_childcareSG.bw, main = "bw.diggle")
plot(kde_childcareSG.ppl, main = "bw.ppl")

11 6.3 Working with Different Kernel Methods

By default, the kernel method used in density.ppp() is gaussian. But there are three other options, namely: Epanechnikov, Quartic and Dics. Let us take a look at what they look like:

par(mfrow=c(2,2))
plot(density(childcareSG_ppp.km,
             sigma=bw.ppl,
             edge=TRUE,
             kernel="gaussian"),
     main="Gaussian")

plot(density(childcareSG_ppp.km,
             sigma=bw.ppl,
             edge=TRUE,
             kernel="epanechnikov"),
     main="Epanechnikov")

plot(density(childcareSG_ppp.km,
             sigma=bw.ppl,
             edge=TRUE,
             kernel="quartic"),
     main="Quartic")

plot(density(childcareSG_ppp.km,
             sigma=bw.ppl,
             edge=TRUE,
             kernel="disc"),
     main="Disc")

Observations: In this dataset, the choice of kernel function has only a minor impact on the overall density plots. The Gaussian, Epanechnikov, and Quartic kernels produce smoother transitions and distribute the density over a broader area. In contrast, the Disc kernel provides a more localized density estimation with sharper boundaries and less smoothness.

12 Fixed and Adaptive KDE

12.1 Computing KDE by using Fixed Bandwidth

Next, we will compute a KDE layer by defining a bandwidth of 600 meter. Notice that in the code block below, the sigma value used is 0.6. This is because the unit of measurement of childcareSG_ppp.km object is in kilometer, hence the 600m is 0.6km.

In this section, we will learn how to derive adaptive kernel density estimation by using density.adaptive() of spatstat.

kde_childcareSG_adaptive <- adaptive.density(childcareSG_ppp.km, method="kernel")
plot(kde_childcareSG_adaptive)

We can compare the fixed and adaptive kernel density estimation outputs by using:

par(mfrow=c(1,2))

plot(kde_childcareSG.bw, main = "Fixed bandwidth")
plot(kde_childcareSG_adaptive, main = "Adaptive bandwidth")

12.2 Converting KDE Output into Grid Object

To achieve the same result, we convert the object to a format suitable for mapping:

gridded_kde_childcareSG_bw <- as.SpatialGridDataFrame.im(kde_childcareSG.bw)
spplot(gridded_kde_childcareSG_bw)

12.2.1 Converting Grided Output into Raster

Next, we will convert the gridded kernel density objects into RasterLayer object by using raster() of raster package.

kde_childcareSG_bw_raster <- raster(kde_childcareSG.bw)

Let us take a look at the properties of kde_childcareSG_bw_raster RasterLayer.

kde_childcareSG_bw_raster
class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.4170614, 0.2647348  (x, y)
extent     : 2.663926, 56.04779, 16.35798, 50.24403  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : layer 
values     : -1.005814e-14, 28.51831  (min, max)

Note that the crs property is NA.

12.2.2 Assigning Projection Systems

To include the CRS information on kde_childcareSG_bw_raster RasterLayer, we will do the following:

projection(kde_childcareSG_bw_raster) <- CRS("+init=EPSG:3414")
kde_childcareSG_bw_raster
class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.4170614, 0.2647348  (x, y)
extent     : 2.663926, 56.04779, 16.35798, 50.24403  (xmin, xmax, ymin, ymax)
crs        : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs 
source     : memory
names      : layer 
values     : -1.005814e-14, 28.51831  (min, max)

Now, the crs property is completed.

12.3 Visualising the output in tmap

Finally, we will display the raster in cartographic quality map using tmap package.

tm_shape(kde_childcareSG_bw_raster) +
  tm_raster("layer", palette = "viridis") +
  tm_layout(legend.position = c("right", "bottom"), frame = FALSE)

values(kde_childcareSG_bw_raster)

Note that the raster values are encoded explicitly onto the raster pixel using the values in “v” field.

12.4 Comparing Spatial Point Patterns using KDE

In this section, we will learn how to compare KDE of childcare at Punggol, Tampines, Chua Chu Kang and Jurong West planning areas.

12.4.1 Extracting Study Area

The code block below will be used to extract the target planning areas.

pg <- mpsz_sf %>%
  filter(PLN_AREA_N == "PUNGGOL")
tm <- mpsz_sf %>%
  filter(PLN_AREA_N == "TAMPINES")
ck <- mpsz_sf %>%
  filter(PLN_AREA_N == "CHOA CHU KANG")
jw <- mpsz_sf %>%
  filter(PLN_AREA_N == "JURONG WEST")

Plotting the target planning areas:

par(mfrow=c(2,2))
plot(pg, main = "Punggol")

plot(tm, main = "Tampines")

plot(ck, main = "Choa Chu Kang")

plot(jw, main = "Jurong West")

12.4.2 Creating owin object

Now, we will convert these sf objects into owin objects that is required by spatstat.

pg_owin = as.owin(pg)
tm_owin = as.owin(tm)
ck_owin = as.owin(ck)
jw_owin = as.owin(jw)

12.4.3 Combining Childcare Points and the Study Area

To extract childcare that is within the specific region for analysis, we can use:

childcare_pg_ppp = childcare_ppp_jit[pg_owin]
childcare_tm_ppp = childcare_ppp_jit[tm_owin]
childcare_ck_ppp = childcare_ppp_jit[ck_owin]
childcare_jw_ppp = childcare_ppp_jit[jw_owin]

Next, rescale.ppp() function is used to transform the unit of measurement from meter to kilometer.

childcare_pg_ppp.km = rescale.ppp(childcare_pg_ppp, 1000, "km")
childcare_tm_ppp.km = rescale.ppp(childcare_tm_ppp, 1000, "km")
childcare_ck_ppp.km = rescale.ppp(childcare_ck_ppp, 1000, "km")
childcare_jw_ppp.km = rescale.ppp(childcare_jw_ppp, 1000, "km")

The code block below is used to plot the four study areas and the locations of the childcare centres.

par(mfrow=c(2,2))

plot(childcare_pg_ppp.km, main="Punggol")
plot(childcare_tm_ppp.km, main="Tampines")
plot(childcare_ck_ppp.km, main="Choa Chu Kang")
plot(childcare_jw_ppp.km, main="Jurong West")

12.4.4 Computing KDE

The code chunk below will be used to compute the KDE of these four planning area. bw.diggle method is used to derive the bandwidth of each area.

par(mfrow=c(2,2))
plot(density(childcare_pg_ppp.km,
             sigma=bw.diggle,
             edge=TRUE,
             kernel="gaussian"),
     main="Punggol")
plot(density(childcare_tm_ppp.km,
             sigma=bw.diggle,
             edge=TRUE,
             kernel="gaussian"),
     main="Tampines")
plot(density(childcare_ck_ppp.km,
             sigma=bw.diggle,
             edge=TRUE,
             kernel="gaussian"),
     main="Choa Chu Kang")
plot(density(childcare_jw_ppp.km,
             sigma=bw.diggle,
             edge=TRUE,
             kernel="gaussian"),
     main="Jurong West")

12.4.5 Computing Fixed Bandwidth KDE

We will set the bandwidth as 250m for comparison purposes.

par(mfrow=c(2,2))

plot(density(childcare_ck_ppp.km,
             sigma=0.25,
             edge=TRUE,
             kernel="gaussian"),
     main="Chou Chu Kang")

plot(density(childcare_jw_ppp.km,
             sigma=0.25,
             edge=TRUE,
             kernel="gaussian"),
     main="Jurong West")

plot(density(childcare_pg_ppp.km,
             sigma=0.25,
             edge=TRUE,
             kernel="gaussian"),
     main="Punggol")

plot(density(childcare_tm_ppp.km,
             sigma=0.25,
             edge=TRUE,
             kernel="gaussian"),
     main="Tampines")

13 Nearest Neighbour Analysis

In this section, we will perform the Clark-Evans test of aggregation for a spatial point pattern by using clarkevans.test() of statspat.

Note

Clarks-Evans Test

The Clark-Evans test is a statistical test used to determine whether points in a given space are randomly distributed, clustered, or regularly spaced.

The Clark-Evans test is a statistical test used to determine whether points in a given space are randomly distributed, clustered, or regularly spaced. It is commonly used in spatial analysis, particularly in nearest neighbor analysis, to understand the pattern of points (such as locations of trees, animals, or buildings) in a study area.

How the Clark-Evans Test Works:

  • The test calculates the average distance from each point in the dataset to its nearest neighbor.
  • This observed average distance is compared to the expected average distance for a random distribution of points in the same area.
  • The result is a Clark-Evans ratio (R):
    • R = 1: Points are randomly distributed.
    • R < 1: Points are clustered together.
    • R > 1: Points are more evenly spaced or regularly distributed.

see Tests for Spatial Randomness - The R Book [Book]

The test hypotheses are:

Ho = The distribution of childcare services are randomly distributed.

H1= The distribution of childcare services are not randomly distributed.

The 95% confident interval will be used.

13.1 Testing Spatial Point Patterns using Clark and Evans Test

clarkevans.test(childcareSG_ppp,
                correction="none",
                clipregion="sg_owin",
                alternative=c("clustered"),
                nsim=99)

    Clark-Evans test
    No edge correction
    Z-test

data:  childcareSG_ppp
R = 0.55631, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)

The test hypotheses are:

\(H_o\) = The distribution of childcare services are randomly distributed.

\(H_1\) = The distribution of childcare services are not randomly distributed.

The 95% confident interval will be used.

13.2 Testing spatial point patterns using Clark and Evans Test

clarkevans.test(childcareSG_ppp,
                correction="none",
                clipregion="sg_owin",
                alternative=c("clustered"),
                nsim=99)

    Clark-Evans test
    No edge correction
    Z-test

data:  childcareSG_ppp
R = 0.55631, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)
Note

Result Interpretation

  • Given that R, the Clark-Evans ratio, is 0.55631, it indicates that the childcare services are clustered rather than randomly distributed. The further R is from 1 (and closer to 0), the stronger the clustering.
  • The p-value is extremely small (much less than the standard significance level of 0.05).
  • At a 95% confidence level, we can reject the null hypothesis of random distribution and accept the alternative hypothesis that the services are clustered.

13.3 Clark and Evans Test: Choa Chu Kang Planning Area

Now, we perform the similar test on individual subzone areas. We will analyse the Choa Chu Kang area first.

clarkevans.test(childcare_ck_ppp,
                correction="none",
                clipregion=NULL,
                alternative=c("two.sided"),
                nsim=999)

    Clark-Evans test
    No edge correction
    Z-test

data:  childcare_ck_ppp
R = 0.95163, p-value = 0.4698
alternative hypothesis: two-sided
Note

Result Interpretation

  • Given that R, the Clark-Evans ratio, is 0.94406, it indicates that the distribution of childcare services is closer to random with minimal clustering. Since R is close to 1 and not significantly less, there is little evidence of clustering or regular spacing.

The p-value is 0.4032, which is much larger than the standard significance level of 0.05.

At a 95% confidence level, we fail to reject the null hypothesis of random distribution. The two-sided alternative hypothesis suggests that the distribution could be either clustered or regularly spaced, but the data does not provide enough evidence to support significant clustering or regularity.

13.4 Clark and Evans Test: Tampines Planning Area

Next, we will analyse the spatial point patterns of childcare centre in Tampines planning area using the Clark-Evans test

clarkevans.test(childcare_tm_ppp,
                correction="none",
                clipregion=NULL,
                alternative=c("two.sided"),
                nsim=999)

    Clark-Evans test
    No edge correction
    Z-test

data:  childcare_tm_ppp
R = 0.79612, p-value = 0.0002337
alternative hypothesis: two-sided
Note

Result Interpretation

Given that R, the Clark-Evans ratio, is 0.79076, it indicates that the distribution of childcare services shows a significant level of clustering. Since R is notably less than 1, there is strong evidence of clustering.

The p-value is 0.0001592, which is much smaller than the standard significance level of 0.05.

At a 95% confidence level, we reject the null hypothesis of random distribution and accept the alternative hypothesis that the distribution is not random. The two-sided alternative hypothesis suggests that the distribution could be either clustered or regularly spaced, but in this case, the significantly lower R value supports the conclusion of clustering.